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We discuss a numerical method to compute the homogeneous solutions of the Teukolsky 
equation which is the basic equation of the black hole perturbation method. We use the 
formalism developed by Mano, Suzuki and Takasugi, in which the homogeneous solutions 
of the radial Teukolsky equation are expressed in terms of two kinds of series of special 
functions, and the formulas for the asymptotic amplitudes are derived explicitly. Although 
the application of this method was previously limited to the analytical evaluation of the 
homogeneous solutions, we find that it is also useful for numerical computation. We also 
find that so-called "renormalized angular momentum parameter", i/, can be found only in 
the limited region of ui for each Z, m if we assume v is real (here, u is the angular frequency, 
and I and m are degree and order of the spin- weighted spheroidal harmonics respectively). 
We also compute the flux of the gravitational waves induced by a compact star in a circular 
orbit on the equatorial plane around a rotating black hole. We find that the relative error 
of the energy flux is about 10"^'' which is much smaller than the one obtained by usual 
numerical integration methods. 

§1. Introduction 

Inspirals of stellar-mass compact objects into a supermassive black hole at galac- 
tic nuclei are expected to be one of the most important sources of the gravitational 
waves for space-based detectors, such as the Laser Interferometer Space Antenna 
(LISA).^) Current best estimate of the number of such event is given by Gair et 
al.^^ They estimated the number of event for the inspirals of WMq black holes into 
10^ Mq suppermassive black holes to be 660 by 3 years observation. By observing 
gravitational waves from such systems, we may be able to obtain information of the 
central black hole's spacetime geometry encoded in multipole moments, and test the 
validity of the no-hair theorem of black hole.^^ We may also obtain astrophysical in- 
formation such like the mass distribution of compact objects in galactic nuclei. The 
optimal detection strategy for such gravitational waves is matched filtering, which 
requires theoretical waveforms to be correlated with the data. Although we may not 
need very accurate waveforms in detection, we will need very accurate theoretical 
waveforms to extract astrophysical information concerning the source. 

To predict the waveforms of extreme mass ratio inspirals, we adopt the black 
hole perturbation approach. In this approach, a compact star is treated as a point 
particle, and the mass of the compact star, fj,, is assumed to be very small compared 
to the mass of the black hole, M, i.e. fi/M <C 1. In this context, the Teukolsky equa- 
tion^^ describes the evolution of a perturbation of the Kerr black hole spacetime. The 
standard approach to solve the Teukolsky equation is based on the Green function 
method. The Green function is constructed using two kinds of homogeneous solu- 
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tions. The solution of the Teukolsky equation is obtained by integrating the Green 
function multipUed by the source term. In the case of a Kerr black hole, the homo- 
geneous solution is calculated usually by the Sasaki-Nakamura equation, which is 
derived with the Sasaki-Nakamura transformation from the Teukolsky equation. The 
Sasaki-Nakamura transformation is a generalization of the Chandrasekhar transfor- 
mation^^ by which we can obtain the Regge- Wheeler equation from the Bardeen- 
Press-Teukolsky equation^)' '''^ in the Schwarzschild case. The Sasaki-Nakamura equa- 
tion is also a powerful formula when we compute the gravitational wave flux induced 
by particles in unbound orbits. 

In the past, there were several works which calculated the flux and waveforms 
of gravitational waves induced by a compact star around a black hole computed 
by the Teukolsky and the Sasaki-Nakamura equations, and its effects to the orbital 
evolution of the star under the influence of radiation reaction. Here, we list some of 
recent works treating bound orbits of the star. Many other works can also be found 
in a review article by Nakamura et al.^) 

Shibata^^ calculated the gravitational waves in the case of circular, equatorial 
orbits around a Kerr black hole. Finn and Thorne^^^ consider the same case, and 
discussed the signal-to-noise ratio of such gravitational waves when detected by 
LISA and the detectability of such sources by LISA. Tanaka et al."*^^) and Cutler 
et al.^^) calculated the gravitational waves in the case of eccentric orbits around a 
Schwarzschild black hole, and discussed orbital evolution of the star. Apostolatos 
et al.^^) consider the slightly eccentric orbits around a Schwarzschild black hole, 
and discussed the stability of the circular orbit under the influence of radiation re- 
action. Kennefick^^) consider the slightly orbits on the equatorial plane around a 
Kerr black hole, and discussed the stability of the circular orbit under the influence 
of radiation reaction. Shibata^^^ calculated the gravitational waves in the case of 
eccentric, equatorial orbits around a Kerr black hole, and discuss the importance 
of the black hole spin and other relativistic effects. Glampedakis and Kennefick^^) 
calculated the gravitational waves in the case of eccentric, equatorial orbits around 
a Kerr black hole, and discuss the orbital evolution of the star under the influence 
of radiation reaction. Shibata^'') calculated the gravitational waves in the case of 
circular, nonequatorial orbits around a Kerr black hole, and discussed the radiation 
reaction effects. Hughes^^)'^^^ calculated the gravitational waves in the case of circu- 
lar, nonequatorial orbits around a Kerr black hole, and discuss the radiation reaction 
effects. 

The source term of the Teukolsky equation is obtained by specifying the orbit 
of the particle. In the case of a Kerr black hole, orbits around it are specified by 
the energy, the z-component of the angular momentum, and the Carter constant. 
In the case of bound orbits, when the orbit is limited to the equatorial plane but is 
not circular, the orbits show exhibit the "zoom- whirl" property as the eccentricity 
becomes larger. When the orbit is no longer limited to the equatorial plane, which 
may be important for sources of LISA, the orbits become more complicated. In such 
cases, we have to trace the orbit for much longer time than the dynamical time of 
the system in order to integrate the source term multiplied by the Green function 
accurately. Further, many side-bands in the spectrum of gravitational waves must 
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be calculated to establish the good accuracy. 

Let us assume for simplicity that LISA will observe gravitational waves at the 
frequency / ~ 10~^Hz for one year. The total cycle of waves is typically Ncyde ~ 
10^. Thus, the relative error for the luminosity AE/E required to establish the 
accuracy for the cycle, Z\A^cycie < 1 is AE/E < 10"^^°^ in the simplest, circular 
orbit cases. For more complicated orbits, the requirement to the accuracy would be 
stronger than this. Although the accuracy, 10~^, is already established in many of 
previous works, it would be very helpful for the future data analysis of LISA if we 
had more efficient and accurate methods to compute the homogeneous solutions of 
the Teukolsky equation to calculate the gravitational wave flux. 

Among approaches to obtain homogeneous solutions of the Teukolsky equation, 
Leaver^^) formulated a method to express them in terms of a series of Coulomb 
wave functions. Although this series is convergent at the spatial infinity, it does not 
converge around the horizon. For this reason, it was suggested that a power series 
expansion around the horizon be used to obtain all of the asymptotic amplitudes 
of the homogeneous solution which are needed in the Green function. Using this 
method, Tagoshi and Nakamura^^-* carried out a high precision computation of a 
homogeneous solution of the Teukolsky equation and obtained high precision nu- 
merical data of the gravitational wave flux induced by a particle in a circular orbit 
around a Schwarzschild black hole. In this calculation, the convergence of the series 
of Coulomb wave functions is very fast, and we can evaluate the series of Coulomb 
wave functions very accurately. Contrastingly, because the convergence of the power 
series expansion is very slow, it is not easy to evaluate it very accurately. In fact, 
the accuracy of the numerical data of the homogeneous solution is limited by the 
accuracy of the power series expansion. 

Later, Mano, Suzuki and Takasugi (MST)^^) formulated a method to express a 
homogeneous solution in a series of hyper geometric functions around the horizon. 
This series has a very convenient property that allows it to be matched with the 
series of Coulomb wave functions around the infinity. Specifically, this property 
is that the three term recurrence relation among the expansion coefficients for the 
hypergeometric case is the same as that in Coulomb case. Owing to this property, it 
is possible to match the two series expansions analytically, and we can express the 
asymptotic amplitude of the homogeneous solution in terms of only the expansion 
coefficients, not in terms of the hypergeometric or Coulomb functions themselves. 
This property is that which distinguishes the MST method from the Leaver method. 
Because we do not need to use a power series expansion to determine the asymptotic 
amplitudes by numerical matching, it is expected that the matching can be done 
very efficiently in the MST method. 

So far, the MST method has been used only in analytic calculations. There is a 
close relation between the Coulomb or hypergeometric series expansion and the low 
frequency expansion of the Teukolsky equation. Here, low frequency usually implies 
small post-Minkowskian or post-Newtonian expansion parameters. We can calculate 
the higher order terms of the post-Newtonian expansion of the gravitational wave 
flux from a black hole^^^'^^^ systematically using this method. Tagoshi, Mano and 
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Takasugi^^) also computed the energy absorbed by a rotating black hole analytically 
assuming that the orbit of the particle is very large. 

In this paper, we use the MST formalism for numerical computation of the 
homogeneous Teukolsky equation. We discuss the numerical method to calculate 
the homogeneous solutions in detail. One of the most important problem involving 
the MST formalism is to determine the so-called " renormalized angular momentum" 
by solving the equation g^v) = 0, which is expressed in terms of continued fractions. 
We investigate the numerical properties of this function g{i^)- We find that, for each 
s,i,m and q {s is the spin index of the Teukolsky equation, i,m are the indices 
of the spin weighted spheroidal harmonics, and q is the Kerr parameter divided by 
the mass of black hole, q = a/M.), there is a maximum value of Mco for which we 
can find v, assuming v to be real. Because the applicability of the MST formalism 
depends on the existence of v, this is a serious problem. Although our preliminary 
investigation suggests that we can find u if we assume it to be complex, more effort 
is needed to establish the numerical accuracy. Therefore, we restrict the region of 
Mcu for which we can find real v. 

Once we have the renormalized angular momentum u, it is straightforward to 
evaluate the expansion coefficients and the series of hypergeometric functions or 
Coulomb wave functions. As expected, we find that the convergence of the series 
of hypergeometric functions and Coulomb wave functions is very fast, and we can 
evaluate the homogeneous solutions very accurately. 

As a test calculation, we compute the gravitational wave flux induced by a 
particle in a circular orbit on the equatorial plane around a Kerr black hole. This is 
because that the numerical data for the energy flux can easily be compared with the 
one obtained by other methods. The accuracy of the numerical data is compared 
with that of previous works. 

This paper is organized as follows. In § 2, we review the Teukolsky formalism 
and the MST method. In § 3, we discuss the numerical method to calculate the 
homogeneous solutions. We present numerical results of the gravitational wave flux 
indTiccd by a particle in circular, equatorial orbits around a Kerr black hole in § 3.5. 
§ 4 is devoted to the summary and discussion. Throughout this paper, we use units 
in which G = c = 1. 

§2. Analytic solutions of the homogeneous Teukolsky equation 

2.1. The homogeneous solutions in series of hypergeometric functions and Coulomb 
wave functions 

The radial Teukolsky equation is given by (see Appendix A) 

^^|:(^^)-nOi^w = rw, (2-1) 

where the potential term V{r) is given by 
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Here A = r"^ - 2Mr + = (r - r+)(r - r_) with r± = M ± y/W^ - a^, K = 
(r^ + a^)uj — ma and A is the eigenvalue of the angular Teukolsky equation. 

In the MST method, the homogeneous solutions of the Teukolsky equation are 
expressed in terms of two kinds of scries of special functions. ^'^^'^^^ One consists 
of series of hypergeometric functions, and the other consists of series of Coulomb 
wave functions. The former is convergent at the horizon and the latter at infinity. 
We match the two kinds of solutions in the overlapping region of convergence. We 
thereby obtain analytical expressions of the asymptotic amplitudes of the solution. 
(See Ref. 25) for a more recent review.) 

First, we present a solution consisting of a series of hypergeometric functions. 
For the incoming solution -RJ^^^^, we define pin by 

Kno. = e-«-(-x)-^-^(^+-)/2(l - xy^'--y^pi^{x). (2-3) 
The function is expressed in a series of hypergeometric functions as 

oo 

Pinix) = anF{n + v + 1 — iT,—n — v — ir;! — s — ie — iT;x), (2-4) 

n=— oo 

where x = ci;(r+ — r)/€K, e = 2Mlo,k = y/l — q^,q = a/M,T = {e — mq)/K and 
F(a, p-jj; x) is the hypergeometric function. 

This expression can be rewritten, using the analytic properties of hypergeometric 
functions, into the form of a series expansion with better convergence properties for 
large |x| as 



where 



= e^^''^(-x)-^-^(^+^)(l - x)^(^+^)+'^ 

r{l-s-ie- iT)r{2n + 2z/ + 1) 



oo 

X 

n=— oo 



r{n + 1^+1- iT)r{n + V + s-ie) 



X (1 - xYF{-n -u-iT,-n-u-s-ie; -2n - 2u; j—^) ■ (^'S) 

Next, we present a solution in the form of series of Coulomb wave functions. Let 
us denote a Teukolsky function by Rc- We define the function fi,{z) through the 
relation 

ifc = .-'-(l-f)-'"-^"'V„(.). (2.7) 
The function fu{z) is expressed in a series of Coulomb wave functions as 

W = E {-i A''^]^'~''K nFn+,{-is - 6, z), (2-8) 

where z = io{r — r-), {a)n = r{a + n)/r{a) and Fj^{r], z) is a Coulomb wave function 
defined by 

Fr,{r,, z) = e-2^z^+i ^^r{2N + 2)^ ^^^ + l-ir^,2N + 2; 2iz). (2-9) 
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Here ^{a,f5;z) is the confluent hypergeometric function, wfiich is regular at z = 
(sec § 13 of Ref. 28)). 

In this method, the expansion coefficients of the scries of hypergeometric func- 
tions and the series of Coulomb wave functions {a„} exhibit the same recurrence 
relation. We find that the expansion coefficients satisfy the three-term recurrence 
relation 

<a„+i + + 7n«n-l = 0, (2-10) 

where 

ieK{n -I- 1/ -I- 1 -I- s -I- ie) (n -I- 1/ -I- 1 -I- s — ie) (n -I- 1/ -I- 1 -I- zr) 



p-^ = -X-s{s + \) + {n + u){n + v + l)+e' + e{e-mq)+'^^^^^^^y 

ieK{n + v — s + ie){n + ly — s — ie){n + v — it) 
^" ~ {n + v){2n + 2v- 1) ' 

We note that the parameter u introduced in the above formulas does not exist in 
the Teukolsky equation. This parameter is introduced so that both series converge 
and actually represent a solution of the Teukolsky equation. We next introduce the 
following quantities: 

L„=^. (2-11) 
We can express and L„ in terms of continued fractions as 

r? — 

rCn. — — 



(3- + a-i?„+i 

7n "n7n+l "n+l7^^+2 



(2-12) 



n+l l~'n+2~ 



< «n-l7n «n-27n-l .2.13) 

The expressions for i?„ and L„ are valid if these continued fractions converge. It 
has been proved^^^ that the continued fraction of the right-hand side of Eq. (2-12) 
converges if and only if the recurrence relations Eq. (2-10) possess a minimal solution 
as n ^ oo. A similar theorem can be proven regarding the converge as n ^ — oo of 
the right-hand side of Eq. (2-13). Because the recurrence relation Eq. (2-10) possess 
minimal solutions as n — >^ ±oo, the continued fractions on the right-hand sides of Eqs. 
(2-12) and (2-13) converge. Although minimal solutions in the limits n ^ oo and 
n ^ — oo do not coincide in general, we can match them by appropriately choosing 
V. Suppose {fn} is ^ solution that is minimal for both n ±oo. It is proved^^-* that 
the following relations are satisfied: 

P — fn _ In "n7n+l '^n+l7n+2 1 4^ 

^n- - av_ Rv _ Ru _ l^'-^^J 

Jn-1 Pn Pn+1 Pn+2 
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This implies the relation 



'y^ 

TV Jn '-'n n—l in in— I in -i r\ 

- -JIT- - lolTZ 'J^ ~ Z • 

Pn Pn-l Pn-2 



RnLn-1 = 1. (2-16) 



If wc choose such that it satisfies the implicit equation for u, Eq. (2-16), for a 
certain n, we can obtain a minimal solution that is valid over the entire range 
— oo < n < oo. For the minimal solution, we have 

lim = hm = (2-17) 

The minimal solution is important for the convergence of the series Eq. (2-4) and 
Eq. (2-8). It can be proved that if we use the minimal solution {/^} for the expansion 
coefficients {a^i}, the scries of hypcrgcomctric functions Eq. (2-4) converges for x in 
the range — oo < x < 0. (In fact, this is true for all complex values of x, except at 
|x| = oo.) It has also been proved that if the expansion coefficients are given by the 
minimal solution, the series Eq. (2-8) converges for z> en or, equivalently, r > r+. 

Instead of Eq. (2-16), we can use an equivalent but practically more convenient 
form of an equation that determines the value of v. Dividing Eq. (2-10) by a„, we 
find 

i^l + + i'J.n-i = 0, (2-18) 

where and are given by the continued fractions Eq. (2-14) and Eq. (2-15) 

respectively. 

2.2. Matching and the asymptotic amplitudes 

Now, we match the two kinds of solutions, Rq and Rq. If we expand solutions 
in powers of5; = l — x = 2:/ {en), we see that both solutions behave like x'^ multiplied 
by a single valued function of x for large | x |, as long as a; > 0. Thus, if we assume 
u > 0, the analytic properties of Rq and Rq are the same. That implies that these 
two solutions are identical up to a constant multiple. We set 

R^ = K^R'S. (2-19) 

Then, by comparing each power of x in the region where both solutions converge, 
i.e. 1 <S| X |< oo, we find 

^ _ e'^'^{2eKy-''-^2-H^r{l - s - ie - iT)r{N + 2u + 2) 
" ~ r{N + U + 1-S + ie)r{N + U + 1 + iT)r{N + v + l + s + ie) 

r{n + N + 2u + I) r{n + u + I + s + ie)r{n + u + l + iT) \ 
' ^ {n-N)\ r{n + i^ + l-s-ie)r{n + iy + l-iT)^''J 

^ {N-n)\{N + 2u + 2)n{i^+l-s + ie)J'' ' ^ ' 
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where can be any integer, and the factor K,^ should be independent of the choice 
of iV. 

Now, we have an alternative expression for the ingoing- wave function R^^^^^, 

Rt^ = K^R^c + if-.-ii?c'"\ (2-21) 

which converges for r > r+, including r = oo. Combining Eq. (2-21) with Eq. (2-4), 
which converges everywhere except at r = oo, we have a complete set of analytic 
solutions for the ingoing-wave function. 

Now, we can obtain analytical expressions for the asymptotic amplitudes Bj^^, 
Bi^^ and Bj^, of R^^ defined in Eq. (A-8). Comparing R}^^ with Eq. (2-4) in 
the limit of r — > r+ and Eq. (2-21) in the limit of r — > oo, we find 



\IjO J 

n=— oo 

ref , ,— 1—23 11/- I A^iT^viy- i av Atlnt 



Rinc , — 1 



where 



B'Cu, = ^~'~''[K. + ie*'^^iC_,_i]A!le"^'^S (2-22) 



r(y + 1 + s — le) ^-^ 



n=— oo 

hoc 



= 2-i-^+-e-^e^^(-+i+^) V (-l)"i^^±i±l^%/-. (2-23) 

n=— oo ^ ' 

We can also obtain an analytical expression for by transforming Rq. We 

note that an analytical property of the confluent hypergeometric function, 

#(a,7,x) = -f^e--'^?Z'(a,7,x) + - a.^i-x). (2-24) 

Fid -01) F{u;) 

Here a = sgn[Im(a;)] and •f'(a, c, x) is the irregular confluent hypergeometric function. 

R'^ = Rl+R'i, (2-25) 



We then rewrite Rq as 



where 



+ Fiv + l + s-ie) ^ ' 



oo 

X 

n=— oo 



E ef!;:^{2zY'F{n + u + l-s + ie,2n + 2v + 2-2iz), (2-26) 



X i''^|^^~'^|" /^(2zWra + I. + 1 + s - ie, 2n + 2z/ + 2; -2i^). (2-27) 
(i/ -I- 1 - s -I- ?e)n 
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For large | a; |, !f'(a,7, x) behaves as 

^{a, P, x) x~" as I a; 1^ oo. (2-28) 

We therefore have 

R'L = A''_z-'^-^'e'^^+'^''^\ (2-29) 

We can see that the functions R'^ and are ingoing-wave and outgoing-wave solu- 
tions at infinity, respectively. In particular, we have the upgoing solution expressed 
in terms of Coulomb wave functions as 

RiL = R-- (2-30) 

Now, we can obtain an analytical expression of the asymptotic amplitude, Cf^^, 
of defined in Eq. (A-9). We find 

/^trans , —l—2sAeln€ av (r) 01 \ 

§3. Numerical methods 

In order to compute the homogeneous solutions of the Teukolsky equation, we 
first compute the eigenvalue A of the spin weighted spheroidal harmonics. This is 
discussed in Appendix B. In this section, we discuss a numerical method to compute 
the homogeneous solution of the radial part of the Teukolsky equation, assuming A 
is given. The eigenvalue A can be computed similar in to v. Although the MST 
formalism was developed for arbitrary values of the spin s, in the rest of paper we 
consider only the case s = —2, because this is important to evaluate gravitational 
waves. 

3.1. Continued fractions 

Because the computation of continued fractions is very important in our numer- 
ical method, we first review Steed's algorithm to compute continued fractions.^''-* Let 
us define a continued fraction by 

^^ = ~^=^o + -r—T—-r—---T-- (3-1) 

Bn bi+ 02+ 03+ bn 



It can be shown that An and Bn satisfy the recurrence relations 

\ ^ f An-1 An-2 \ ( \ / A_i \ / 60 1 

Bn / I Bn-l Bn-2 / I a„ ' [ Bq B^i I 1 



(3-2) 



Here, we introduce the quantity Dn = Bn-i/Bn- Then, the above recurrence relation 
of Bn is identical to 
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We can rewrite the right-hand side of Eq. (3-1) in terms of D„ as 

^ bnAn—l + flfiVln— 2 

Bn bnBn-1 + anBn-2 ' 

_ bn{An,-i/ Bn-i) + an{An-2Bn-2/ Bn-2Bn-l) 

bn + 0"n{Bn-2/Bn-l) 
= {hn-lbn + hn-2Dn-ian)Dn. (3-4) 

Furthermore, from Eq. (3-3), we have the relation anDn-iD^ = 1 — 6„D„, which 
then yields the difference between hn and hn-i-, 

Ahn = hn — hn-l, 

= {bnDn - l)(/ln-l - hn-2), 

= {bnDn-l)Ahn-l. (3-5) 

Steed's algorithm is summarized as follows. We first set ho = bo, Di = 1/bi, Ahi = 
ai/61 and hi = ho+Ahi. We then compute hn beginning with n = 2 until | Ahn/hn \ 
is smaller than the required accuracy. 

3.2. Determination of the renormalized angular momentum v 

The parameter v is determined as a solution of Eq. (2-18). In particular, we set 
n = and obtain 

g{v) = i5l + a^i?i + 7o"L_i = 0, (3-6) 

where R\ and L_i are expressed by the continued fractions Eqs. (2-14) and (2-15), 
respectively. 
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Fig. 1. g{v) forl = 2,m = 2. 



In order to search for a root of the implicit equation g{i>) = 0, we can use various 
numerical techniques. Among them, we adopt Brent's algorithm (e.g. Ref. 31)). In 
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Fig. 3. g{v) for £ = 4, m = 4. 

order to search for a root of implicit equations, we need an initial value of the search 
which is not very far away from the root of the equation. In the case that Ma; 
is small, we can use an analytic expression of v as the initial value. Among the 
infinite number of roots, there is an analytic expression of v in the form of a series 
of e = 2Mlo given by 



1 



2£ + l 



+ 



+ 1) 



„2l2 



{u - i)2e{2e + 1) 



£{£ + 1) (2£ + l)(2£ + 2)(2^ + 3) 
e2 + 0(e3). 



(3.7) 
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Fig. 4. g{i') for £ = 5, m = 5. 

When Muj is not very large (less than 0.36 in the case q = 0), this expression 
is useful as an initial value for the root search algorithm. Note that the analytic 
expression for u truncated at O(e^) is always slightly less than £. Thus, in practice, 
it is sufficient to search for a root in the range £ — 1/2 < u < £. We can see in Figs. 
1-4 that the function g^i') is smooth in this range. Although numerically we find 
poles at half integer values of v, it is not difficult to find a root in this region. In 
fact, we can find a root as long as Moo < 0.36 in the case £ = m = 2 and q = 0. 

However, the situation is different when Mlv becomes larger. In this case, u 
approaches i — 1/2, and beyond a certain value of Muj, it becomes impossible to find 
roots v. The maximum values of Mui for which real v can be found in the region 
£ — 1/2 < V < £ are listed in Table I. The maximum values depend on £,m and 
q. This maximum value of Mco increases with £. The maximum value is larger for 
q > than for q < 0. The behavior of u beyond Mcu in Table I depends on £ and 
m. In some cases, we can find as a real value. However, in other cases, there is no 
root of Eq. (3-6) for any real value of u. 

We note that we have attempted to find a complex root u. Our preliminary 
investigation suggests that there is such a root in the complex value of when we 
can not find a real root. However, more work is needed to confirm this, and we leave 
this as a future project. In the rest of this paper, we concentrate on the case of real 
v and the region of e in which we can find such a root. The values of v for various 
values of Muj = e/2, £, and q arc plotted in Fig. 5. 

Note that there are some roots at integer and half integer values of f that are 
independent oi £,m,q, and uj. Although they are solutions of = 0, they do not 
connect the two minimal solutions in the limits n ^ ±00. 

When u is a root of Eq. (3-6), the quantities v + k {k = ±1, ±2, • • • ) are also 
roots of Eq. (3-6), since u only appears as + n in the continued fractions of Eq. 
(3-6). This fact can be checked numerically for k = ±1. However, as we can see 
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Fig. 5. I' as a function of Muj for £ = m = 2 to £ = m = 5 for g = 0, 0.9 and —0.9. 
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Table I. Maximum values of Muj for which real f is found. 



£ m 
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q = -0.9 
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in Figs. (l)-(4), the function g{i') around the roots ± 1 is a very steep function. 
When \k\ is larger than 2, the slope of g{u) becomes much steeper, and it becomes 
difficult to find the root u + k. 

3.3. Expansion coefficients 

Once we have v, it is straightforward to evaluate the expansion coefficients using 
the three-term recurrence relations Eq. (2-10). As discussed in § 2, the expansion 
coefficients {a^} must be the minimal solution {/n}- It is well known that the 
minimal solution cannot be computed numerically with forward recursion of the 
three-term recurrence relations from n = to ±cx). This is because a small numerical 
errors in the minimal solution /„ contain dominant parts of the solution of the 
recurrence relation. For this reason, the numerical solution obtained using forward 
recursion will be dominated by errors after several recursions, since the dominant 
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Table 11. ;/ tor AI^ = O.i. 
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solutions grow very rapidly. 

In such a situation, we can use the continued fractions given in Eqs. (2-14) 
and (2-15) to evaluate the minimal solution. First, wc set the initial value of as 
/o = 1. Wc evaluate fi from /i//o, which is evaluated by the continued fraction, the 
right-hand side of Eq. (2-14) with n = 1. We then obtain /i as /i = Rifo- In the 
same way, we can evaluate /n (n > 1) from fn-i and Rn recursively. For n < 0, we 
use the same algorithm and evaluate /„ from fn+i and L„. 

3.4. Homogeneous solutions 

Given the eigenvalue of the spin weighted spheroidal harmonics A, the rcnormal- 
ized angular momentum zv, and the expansion coefficients /„, it is straightforward to 
compute the asymptotic amplitude of the homogeneous solutions, Eqs. (2-22) and 
(2-31). It is also straightforward to compute the homogeneous functions using Eq. 
(2-3) or Eq. (2-30). We found that the convergence of the series of hypergeometric 
functions or Coulomb wave functions and that of the formulas for the asymptotic 
amplitude are very rapid. This is because the expansion coefficients {/^}, which 
constitute a minimal solution of the three-term recurrence relation, decrease very 
rapidly as |n| becomes large, as can be seen from Eq. (2-17). For example, we only 
need to carry out the summation in Eq. (2-4) from n = to n '-^ ±20 to have an 
accuracy 10^^^ in the case e < 0.4. 

The hypergeometric function F{a, /3; 7, x) is computed using the transformation 
formula 

r(y)r(3 - a) 1 
F{a, (3; 7, x) = {1 - x)'" ^ '-F{a, 7 - /?; « - /? + 1; ) 

- -'-' rwrTffl ^'^'^- 1; T^). (3-8) 

and using the Gauss hypergeometric series for the hypergeometric functions on the 
right-hand side of Eq. (3-8). For the gamma function, we use the routines available 
on the world wide web.^^^ 

3.5. Gravitational wave luminosity 

In order to check the accuracy of our numerical code, we calculated the grav- 
itational wave flux from a point particle in circular orbits on the equatorial plane 
around a Kerr black hole. The formula for the luminosity is given in Appendix A. 

Our computations presented here were done using a double precision code. In 
Tables III-VI, we hst the luminosities for ro = 6, 10, 100 and 1000 from ^ = 2 to 6 
and for g = and q = ±0.9. We also list the total luminosity for various values q 
and ro in Table VII. 

We compared these results with those of Tagoshi and Nakamura^^^ in the Schwarzschild 
case q = 0. The results are given in Table VIII for cases ro = lOM. We find that 
our results agree with them with relative error 10~^^ - 10~^^. Because the estimated 
accuracy of Ref. 22) is more than 20 significant figures, we estimate that the accu- 
racy of our code is about 13-14 significant figures. We also compare our results with 
those given by Kennefick^^^ in the Kerr case q ^ 0. The results are given in Table IX 



New Numerical Methods to Evaluate the Homogeneous Teukolsky Equation 17 

for cases rg = 3M and q = 0.998. We find that our results agree with them with rel- 
ative error 10~^ - 10~^. This is consistent with the estimated accuracy in his works. 
Although we do not have a quantitative estimate of the accuracy for the luminosity 
in the Kerr case, because u is computed with accuracy in the Schwarzschild case, we 
expect that the accuracy of our results is also about 13-14 significant figures in the 
Kerr case. 



Table III. The gravitational wave luminosity for ro = 6M. 
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-6 


6.6947435866 
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3 


3.4798321412 
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1.4534938751 
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-4 


8.0343009373 
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-5 


3 
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-8 
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-5 
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-8 
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X 
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-8 


1.7547382460 
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10" 


-8 


6.1994559709 


X 


10" 


-9 


4 


1 


4.5312182607 


X 
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11 


7.0534448975 


X 
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12 


9.3343722027 


X 
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13 


5 


5 
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2.5069319096 
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19 



§4. Summary and discussion 

In this paper, we described numerical methods to compute the homogeneous 
solutions of the Teukolsky equation. We used the MST formalism, in which the 
homogeneous solutions of the Teukolsky equation are expressed in terms of series of 
hypergeometric functions and Coulomb wave functions. We found that the renor- 
malized angular momentum v can be found only for a limited values of Mlu. When 
Mu becomes large, u approaches £—1/2, and over a certain range of values of Mu, 
which depends on i, m and q, we could not find real v. Our preliminary investiga- 
tions suggest that becomes complex when Muj becomes sufficiently large. However, 
because further investigation is needed to confirm this, we continue to work on this 
point. 

In the region of Muj in which we can find real ly, we found that the convergence of 
the series of hypergeometric functions or Coulomb wave functions is very rapid. This 
is because the expansion coefficients {/^}, which constitute a minimal solution of the 
three-term recurrence relation, decrease very rapidly as |n| becomes large. Thus, we 
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Table IV. The gravitational wave luminosity for ro = lOM. 
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X 


1 0' 
lU 


-6 


5.0413559470 


X 


lU 


-6 


Q 
O 


9 


1.2628541265 


X 




-7 


4;. / yoy J-U4:U J-U 
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1.1385202800 


X 


10" 


12 


3.8193532372 


X 


10" 


13 
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concluded that the MST formahsm is a powerful method to compute the homoge- 
neous solutions numerically. By comparing the numerical data for the gravitational 
wave luminosity emitted to infinity in the case that a particle display circular orbits 
in the equatorial plane with the results of previous works, we estimated the accuracy 
of our code to be about 13-14 significant figures in the double precision computation. 
This accuracy will be sufficient as a Green function of the Teukolsky equation in the 
computation of the templates for the data analysis of LISA. Currently, the accuracy 
is limited not by the MST formalism itself but by the accuracy of the evaluation of 
the hypergeometric functions. 

In the near future, we will investigate the question of whether complex v exist 
in order to broaden the region of Ma;, in which this method can be applied. We 
also hope to improve the accuracy of the evaluation of the hypergeometric func- 
tions and Coulomb wave functions and to apply this method to the computation of 
gravitational waves from a compact star in the case of more generic orbits around a 
supermassive black hole. 
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Table V. The gravitational wave luminosity for ro = lOOAf . 
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Appendix A 

Teukolsky Formalism, 

A.l. Teukolsky equation 

In terms of the Boyer-Lindquist coordinates (t, r, 0, (f)), the metric of a Kerr black 
hole is expressed as 

ds^ = -^{dt-a sin^ edcPf + [(r2 + a^)^ -adt]^ 

+ ^dr^ + Ede^ = g^pdx'^dx'^. (A-1) 

where E = r"^ + cos^ 9 and A = — 2Mr + a^. In the Teukolsky formalism,^) 
the gravitational perturbations of a Kerr black hole are described by the Newman- 
Penrose quantity 1^4 = —Ca/3-ysn"rn^n'^rn^, where Ca/3j5 is the Weyl tensor and 

m° = [msin e, 0, 1, V sin 6] /{V2{r + ia sin 9)), 

rrf = [-msin6',0,l,-z/sin6']/(V2(r-msin6')). (A-2) 

The Weyl curvature component '04 contains all the information regarding the 
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Table VI. The gravitational wave luminosity for ro = lOOOM. 
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gravitational radiation. Teukolsky showed that if we carry out Fourier-harmonic 
decomposition of p'^^tp^, we can separate the Teukolsky equation into a radial part 
and an angular part, or 

£,Tn 

AttET = Y,j cL;e--*+^-'^ ^2SZ{0)Tim.{r) , (A-3) 

where p = (r — ia cos 9)~^. 

The radial function Rimui{i") and the angular function _25|^(^) satisfy the fol- 
lowing Teukolsky equations: 



1 


d 










sint 





(m — 2 cos ( 



^2 



+Aa(jo cos — 2-1- 2mau + A 

The potential is given by 



sin^ I 



(A-5) 



V{r) = ^ ^ h Stojr + A, (A-6) 
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Table VII. The gravitational wave luminosity, up through i = 6, for various q and orbital radii 
ro/M. 
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where K = [r^ +a'^)u} — ma and A is the eigenvalue of -2S^{6). The angular function 
-2S'^{9) is called a spin-weighted spheroidal harmonic with spin weight —2. It is 
usually normalized as 



p7\ 

Jo 



.2SZ\'^smede = 1. (A-7) 



In the case of a Kerr black hole, analytic form of the spin-weighted spheroidal 
harmonics and the eigenvalue A are not known. However, in the case of a Schwarzschild 
black hole, the spin-weighted spheroidal harmonics reduce to the spin-weighted 
spherical harmonics, and its eigenvalue is A = £{i + 1) — s{s + 1).^^^ 

We solve the radial Teukolsky equation using the Green function method. For 
this purpose, we introduce two kinds of the homogeneous solutions of the radial 
Teukolsky equation. 

Din I ^imw ^ ^ lor r ^ r+, 

^Imw y ^3^rcf^gi<^r* + ^-l^inc^g-it^r* for r ^ +00, ^ ' 

( r'^P pikr* _|_ /\2/^ref -ikr* f„ „ 
„up I l^irruj(^ + ^ '-'Imoj^ tOT r ^ r+, .. . 

^irrn^ \ r^Cf^^^e^''* for +00, ^ ' 
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-9 


1.4921162748527967 


X 


10" 


-9 


1 


91 


X 


10" 


14 


5 


5 


1.5241547645798743 


X 


10" 


-7 


1.5241547645799033 


X 


10" 


-7 


1 


90 


X 


10" 


14 


6 


1 


3.5977953599117842 


X 


10" 


20 


3.5977953599117314 


X 


10" 


20 


1 


47 


X 


10" 


14 


6 


2 


1.9763689535200328 


X 


10" 


16 


1.9763689535200461 


X 


10" 


16 


6 


73 


X 


10" 


15 


6 


3 


2.1238827476369010 


X 


10" 


13 


2.1238827476368560 


X 


10" 


13 


2 


12 


X 


10" 


14 


6 


4 


4.6633398847411170 


X 


10" 


11 


4.6633398847412050 


X 


10" 


11 


1 


89 


X 


10" 


14 


6 


5 


2.4746386947271652 


X 


10" 


10 


2.4746386947272387 


X 


10" 


10 


2 


97 


X 


10" 


14 


6 


6 


2.5182131568101545 


X 


10" 


-8 


2.5182131568101642 


X 


10" 


-8 


3 


82 


X 


10" 


15 


7 


1 


3.2913629491589289 


X 


10" 


23 


3.2913629491588661 


X 


10" 


23 


1 


91 


X 


10" 


14 


7 


2 


9.0841508908488114 


X 


10" 


19 


9.0841508908487538 


X 


10" 


19 


6 


34 


X 


10" 


15 


7 


3 


2.0373627509685863 


X 


10" 


16 


2.0373627509686016 


X 


10" 


16 


7 


53 


X 


10" 


15 


7 


4 


6.9940936502071495 


X 


10" 


14 


6.9940936502074081 


X 


10" 


14 


3 


70 


X 


10" 


14 


7 


5 


1.0340989127935083 


X 


10" 


11 


1.0340989127934911 


X 


10" 


11 


1 


66 


X 


10" 


14 


7 


6 


4.0679948091711724 


X 


10" 


11 


4.0679948091710895 


X 


10" 


11 


2 


04 


X 


10" 


14 


7 


7 


4.2345226712846746 


X 


10" 


-9 


4.2345226712847762 


X 


10" 


-9 


2 


40 


X 


10" 


14 



where k = lo — ma/2Mr^ and r* is the tortoise coordinate defined by 



f dr* , 

= J d^'''-' 
-I 



dr, 



2Mr+ r — r_i_ 2Af r_ r — r_ , . , 

= r + ^In — — t hi — — . (A-10) 

r+ - r_ 2M r+ - r_ 2M ^ ' 

Here r± = M ± - a"^. 

Then, a solution of the radial Teukolsky equation that is ingoing at the horizon 
and outgoing at the infinity can be written as 



,(-•) / 

Jr 



Z\2(r') 
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Tabic IX. Relative error of the energy flux between Kennefick's results^^' and our results for ro = 

:^AL q = 0.998. 
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1 


Kcnncfick 




This paper 








Absolute values 


'ji relative 


2 


1 


6.328361227 x 10" 


-6 


6.3283737145864869 


X 


10" 


-ti 


1.97 


X 


10" 


b 


2 


2 


7.335124682 x 10" 


_3 


7.3350910900390372 


X 


10" 


_3 


4.58 


X 


10" 


6 


3 


1 


3.451784675 x 10" 


-8 


3.4517924227877552 


X 


10" 


-8 


2.24 


X 


10" 


6 


3 


2 


5.645895586 x 10" 


-(i 


5.6458902861747337 


X 


10" 


-6 


9.39 


X 


10" 


7 


3 


3 


2.208481088 x 10" 


-3 


2.2084763695700243 


X 


10" 


-3 


2.14 


X 


10" 


6 


4 


1 


8.842632418 x 10" 
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8. 842634 1846420484 


X 


10" 


12 
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X 


10" 


7 


4 


2 
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9.6507761069096606 
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10" 
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10" 


6 


4 


3 
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3. 1572084515052529 
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10" 
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X 


10" 


6 


4 


4 
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10" 
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X 


10" 


6 


5 


1 


2.226601380 x 10" 


14 


2 . 226600358731 7006 


X 


10" 


14 


4.59 


X 


10" 


7 


5 


2 


7.623615371 x 10" 


11 


7.6236253541558626 


X 


10" 


11 


1.31 


X 


10" 


6 


5 


3 


9.637713992 x 10" 


-8 


9.6377094209411545 


X 


10" 


-8 


4.74 


X 


10" 


7 


5 


4 


1.524036288 x 10" 


-6 


1.5240423964556889 


X 


10" 


-6 


4.01 


X 


10" 


6 


5 


5 


3.044015502 x 10" 


-4 


3.0439980326075501 


X 


10" 


-4 


5.74 


X 


10" 


6 


6 


1 


3.850986300 x 10" 


18 


3.8509863211419323 


X 


10" 


18 


5.49 


X 


10" 


9 


6 


2 


5.736053080 x 10" 


13 


5.7360410503468190 


X 


10" 


13 


2.10 


X 


10" 


6 


6 


3 


1.315574802 x 10" 


10 


1.3155762424734953 


X 


10" 


10 


1.09 


X 


10" 


6 


6 


4 


6.788432814 x 10" 


-8 


6.7884404091806854 


X 


10" 


-8 


1.12 


X 


10" 


6 


6 


5 


6.900330654 x 10" 


-7 


6.9003809076206216 


X 


10" 


-7 


7.28 


X 


10" 


6 


6 


6 


1.219262418 x 10" 


-4 


1.2192534271212692 


X 


10" 


-4 


7.37 


X 


10" 


6 


7 


1 


5.063200000 x 10" 


21 


5.0631852897823082 


X 


10" 


21 


2.91 


X 


10" 


6 


7 


2 


3.106010424 X 10" 


16 


3.1060171837611518 


X 


10" 


16 


2.18 


X 


10" 


6 


7 


3 


1.769856886 x 10" 


12 


1.7698525192439268 


X 


10" 


12 


2.47 


X 


10" 


6 


7 


4 


1.284621326 x 10" 


10 


1.2846205001945888 


X 


10" 


10 


6.43 


X 


10" 


7 


7 


5 


4.051197224 x 10" 


-8 


4.0511763345368075 


X 


10" 


-8 


5.16 


X 


10" 


6 


7 


6 


3.023851152 x 10" 


-7 


3.0238471412501713 


X 


10" 


-7 


1.33 


X 


10" 


6 


7 


7 


4.993950576 x 10" 


-5 


4.9939392339609413 


X 


10" 


-6 


2.27 


X 


10" 


6 



where Wimu is the Wronskian given by 

W,^ = W[A-^'^R^, ^-"^RZJ = '^^^CfZ'BZ^- (A-11) 
The asymptotic behavior at the horizon is given by 



^zl^A\-^^^\ (A-12) 
The asymptotic behavior at infinity is given by 

In this paper, we focus on the gravitational wave fiux from a point particle in 
circular, equatorial orbits around a Kerr black hole. In this case, ^/mo; in Eq. (A- 13) 



24 



R. Fujita and H. Tagoshi 



takes the form 

n 

where ujn = mM^/^/{r^/^ + aM^/^). 

Then, the time-averaged flux (luminosity) radiated to infinity is given by 

^ = y ' ' ■ (A-15) 

l,m,n 

A. 2. Source term 

The source term of the Teukolsky equation, Eq. (A-4), is given by 

/oau) 
df2dtp-'>p-HB', + S-)e-mc,+^a,tzH^, (A.16) 

where B2 and B| are given by 

B'^ = ~p^pL_i[p-^Lo{p-^r'Tnn)] 

^ p^-pA^L^^ip-^p'j+ip-^p-^A-^Trnn)], 



2V2 

B2 = —p^'pA^J+ip'^J+ip-^pTrr.rn)] 

- ^p'pA'j+[p-Y^-'L^i{p-'p-'T^n)] , (A-17) 

TTl 

jCs = de -\ — awsin^ + scot0, (A- 18) 

sine* 

J+ = dr + iK/A ; K = {r'^ + a'^)u - ma, 

and r„„, T-jfin, Tjf^ are the tetrad components of the energy momentum tensor 
{Tnn = Tfj_yn^^n'^ etc.). Here the bar denotes the complex conjugation. 

We consider T^^ of a monopole particle of mass ji. In this case, the energy 
momentum tensor takes the form 



Ti''' =pL j dT 



dzf'dz" 5^^\x-z{t)) 
dr dr \f—9 



where = dz^ /dr and z^ = (t, r{t), 0{t), (p{t)) is a geodesic trajectory and r = r(t) 
is the proper time along this geodesic. The geodesic equations in the Kerr geometry 
are given by 

2 2 

jj^ = [E{r'' + a^) - aL,] - a [aE sin^ 9 - L,] , 
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= ±Vr, 

dr 

dr 



where E and are the energy and the z-component of the angular momentum of a 
test particle, respectively. Also, we have 

R = [{r^ + a^)E - aL,] ^ - A[r^ + (L, - aEf + C], 



= C- 



(1 - + ^ 

sm^ 



cos^ 6, 



where C is the Carter constant of a test particle. Note that E, and C are measured 
in units of /i. If they are expressed in the standard units, E, and C in the above 
equations must be replaced by E/ jjl, L^/ jj, and C j [p- . 

Using Eq. (A- 19), we obtain the tetrad components of the energy momentum 
tensor as 



T — 



sm( 



Mr-r{t))S{Q-Q{tm4>-m\ 



(A-20) 



Trnn = f^^Sir - r(t)),5(^ - 9{t))6icj> - cf>{t)), 
sm a 



(A-21) 



fi^sir - ritme - e{t))d{cj> - m , 

smt/ 



where 

Cnn 



1 



E{r^ + a^) -aL^ + E 



dr 



n 2 



P 



2Eu^ 



2yt 



isinOiaE ^ + 

sin^fc'/ dr 



dr 
d^ 



sm 



(A-22) 



i sin 6 ( aE 1— ^ + E 



dT_ 
(A-23) 



Substituting Eq. (A- 17) into Eq. (A- 16) and performing integration by parts, we 
obtain 



Ten 



dt / d^e*'^*-*™^(*) 



\/27r J-co Jo 
+ ^ {4SZ + ia{p - p) sin 9 SZ) 
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xJ+{Cfnnp-^p-^A-H{r - r{t))S{e - e{t))} 
+ ^Ll{p^SZ{p'p-'),r}CfnnAp-'p-'Sir - r{t))6{e - e{t)) 

- ]p^A^SZJ+{P''MPP~'Cfnfn6ir - r{t))5ie - e{t)))}] , (A-24) 



where 



Cl = da- 



rn 



+ auj smO + s cot I 



smf 



and SZ denotes _2>S'^^(^) for simplicity. 

The integration over can be performed directly. It yields 



Hmuj — P' 



where 



(^nnO + ^mnO + Ammo)S{r - r{t)) 



+ {{Arnnl + ^mml)<^(r - r{t))} ^ + {A^m2'5(r - r{t))} 



A, 



nnO 



Cnnp-'p-'Lt{p-'Ltip''SZ)h 



(A-25) 



9=9{t) 



, (A-26) 



(A-27) 



Crnnp-'[{LtSZ) + P + p) 



-asm9{t)SZ§{p-p) 

K\ 



=p-'-pC^^SZ[-^Q ^^~^ + 2^p^ 
p-^Crnn[LtSZ + iasiu9{t){-p - p)SZl 



K 



2 



p-^pCrnrnSZ{^^+ P) 



27r 

1 



P pC-mmSZ- 



Inserting Eq. (A-26) into Eqs. (A- 12) and (A- 13), we obtain as 

/oo 
-oo 



? trans 



II Rt 



(A-28) 
(A-29) 
(A-30) 

(A-31) 
(A-32) 

(A-33) 



where 



rvrX) 



p 



(A-34) 



■'Emijj 



■^^ma;{^""0 + ^mnO + ^mmo} 
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(A-35) 



dr 



dr'^ 



r=r(t),e=e{t) 



■JOO 



-^^mu){^"wO + -^mnO + ^mmo} 
J pin j2 rjin 



dr 



dr^ 



mm2 



r=r{t),e=e{t) 



(A-36) 



Appendix B 

The Spin- Weighted Spheroidal Harmonics 



In this section, we review the formahsm to represent the spin-weighted spheroidal 

harmonics in a scries of Jacobi polynomials based on Ref. 35). We also discuss a 
method to implement this formalism in the numerical computation. 
We first transform the angular Teukolsky equation as 



m + s + 2msx 



2six + E 



0, 



(B-l) 



where ^ = a x = cos and E = \ -\- s{s-\-l) — a^up- + 2amL0. 

The angular function sSi^{x) is called the spin- weighted spheroidal harmonics. 
Equation (B-l) is a Sturm-Liouville type eigenvalue equation with the boundary 
conditions of the regularity at x = ±1. The eigenvalue E depends on i for fixed 
parameters s, m and au. Therefore we express the eigenvalue as sE'P{^). When au = 
0, sSi^{x) is reduced to the spin- weighted spherical harmonics, and the eigenvalue 
sEJ^i^) becomes ^(^ + 1).34) 

The differential equation Eq. (B-l) has singularities at x 
transform the angular function as 



±1 and oo. We 



1 



1 + x 



where a 



2 J \ 2 

|m + s| and (3 = \m — s\. Then, Eq. (B-l) becomes 
(1 - x^) sUlUx) + [/?-« - (2 + a + P)x] sUUx) 



(B-2) 



+ 



2 \ 2 ' ' ' 
C [-2(1 - x^) sU^x) + {a + (3 + 2s + 2)x sUim{x) 
-{(,+l3-a)sUim{x)]. 



(B-3) 



When ^ = 0, the right-hand side of Eq. (B-3) is zero, while the left-hand side 
becomes the differential equation satisfied by the Jacobi polynomials. 



(1 - x^) P^'^^^y'ix) + [/?-« - (a + /? + 2)x] pj^^^^yix) 

+r(r + a + P + l) P^'^^^^ {x) = 0. 



(B-4) 
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In this limit, the eigenvalue in the equation (B-3) becomes £{i + 1), where 

£ = r + {a + (3) /2 = r + max(| m |, | s |). 

In this paper, we use the Jacobi polynomials defined using the Rodrigue's formula 

by 



\0+n 



fB-5) 



Now, we expand sUim{x) in a series of Jacobi polynomials: 

oo 
n=0 

The expansion coefficients s-^j^ (0 satisfy the recurrence relations 

«^°^4i^(0 + /?^°^A^S(e) = 0, (B.7) 
sAt"\0 + + 7^"^ sAt'\0 = 0, (n > 1), (B-8) 

4C(/) + n + l){n + 3 + l){n + (n + J)/2 + 1 - s) 



where 



(2n + a + /3 + 2)(2n + a + /3 + 3) 

/?<"' = ,Er«) + «^ - (" + ^) (" + ^ + 1) 

2^s(a-/3)(a + /3) 



(2n + a + /?)(2n + a + /3 + 2) ■ 



(n) ^ 4Cn(n + a + ^)(n + (a + /3)/2 + 5) 

' (2n + a + /3- l)(2n + a + /3) ' ^ ' 

The method to determine the eigenvalue gEf^ is similar to that in the case 
of renormalized angular momentum presented in § 3.2. The three-term recurrence 
relation Eq. (B-8) has two independent solutions, which behave for large n as 

An) const. (-0" 

^(^) ~ r(n + (a + /3 + 3)/2-.)' ^^-''^ 
^[2) ~ const.rr(n + (a + /3 + l)/2 + s). (B-11) 

The first one, ^[ij* > is the minimal solution, and the second one, ^pj' ^ dominant 
solution, since lim„_>oo ^(1) M(2) ~ ^- -^"-"^ ^'^^ ^(2) ' ^'^^^^ coefficients increase 

with and the series Eq. (B-6) diverges for all values of x. In the case of this 

series converges. Thus, we have to choose ^^^j* in the series expansion Eq. (B-6). 
Next, we introduce the following: 

AT' 

Pri = 77^, Ln = — qiY- (B-12) 
^(1) ^(1) 
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The ratio Rn can be expressed as a continued fraction: 

Q,(n)^(n+l) Q,(n+l)^(n+2) 

We can also express L„ as a continued fraction: 



/5{n)_ /3(n-2)_ ^(1)_ ^(0) " (^'l^) 

This expression for i?„ is vahd if the continued fraction converges. By using the 
properties of the three-term recurrence relations (Ref. 29), p. 35), it is proved that 
the continued fractions Eq. (B-13) converge as long as the eigenvalue gEl^ is finite. 
Now, we divide Eq. (B-8) by the expansion coefficients s^l^ ^^'^ obtain 

hi sE^) = + a^'^^Rn+i + 7^"^i^n-i = 0. (B-15) 



We replace Rn+i and Ln-i by the continued fractions given in Eqs. (B-13) and 
(B-14). We can determine the eigenvalue gEJ^ as a root of the equation Eq. (B-15), 
where I = n + [a + (3) /2. Because the equations in Eq. (B-15) with different values 
of n are independent, they give different values of gEJ^. We thus obtain each gEJ^ 
for different (. even if the recurrence relation Eq. (B-8) does not contain I explicitly. 

As in § 3.2, we adopt Brent's algorithm^-*^) in order to determine s-E^. In the 
case that auj is not large, we can use an analytical expression of gE'^ as the initial 
value of the root search algorithm. An analytical expression of gE'^ in a series of 
powers of ^ = ao; is given by 

sEf = £{£ + 1) - + [H{e + 1) - H{e) -l]e + Oif), (B-16) 

where 

^^^^ - (2£-l)£3(2£ + l) • (^-^^^ 

The function h{ sE'/'') is usually monotonic, and therefore it is very easy to find the 
root of Eq. (B-15). 

After we obtain the eigenvalue sEY^, we can determine all the coefficients. The 
coefficient for n = rig = I — (q + /3)/2 is usually the largest term. The ratio of other 
terms to those for n = n^, i.e. A^^^^/A^^^\ can be determined using Eqs. (B-13) and 
(B-14) from n = to n = or n = oo. 

The factor A^^^^ is determined by the normalization condition. For this, we 
introduce a new function sVimix) through 

sSZix) = e-«- (^—^\ " " sVUx). (B-18) 
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From Eqs. (B-2) and (B-18), we find 

sVlm{x) = exp(2^x) sUlm{x). 



(B-19) 



We insert Eq. (B-18) into Eq. (B-1) and find that sVlm{x) satisfies the differential 
equation 

(1 - x^) ,V^'4(x) + [/3-a-(2 + a + P)x] ,V^'„(x) 



+ 



a + P f a + P 



+ 1 



^ [2(1 - .t2) ,VUx) - (a + /? - 2s + 2)x sVi^{x) 
~{^-P + a)sVimix)]. 



(B-20) 



In the same way as in the case of sUim{x), we expand s^m(3j) i^i ^ series of 
Jacobi polynomials: 



(B-21) 



n=0 



The expansion coefficients sB^i^{i) satisfy the recurrence relations 

"^°^-<^(0+/3(°^i?S(0 = o, 

sBt^'^ (0 + /3(") (e) + 7^") sBt'^ (0 =0, (n > 1) 



(B-22) 



where 



a 



(n) = _ 



4^(n + g + l)(n + /? + l)(n + (a + /3)/2 + 1 + s) 
(2n + a + /? + 2)(2n + a + ^ + 3) 



H ^ 1- 1 



2^s(a-/3)(a + /3) 



(2n + a + /3)(2n + a + /? + 2) 
4^rz(n + g + /3)(n + (g + /3)/2 - s) 
{2n + a + f3-l){2n + a + (3) 



(B-23) 



In order for the series Eq. (B-21) to converge, the coefficients sBi'^ must con- 

(n) 

stitute the minimal solution of the recurrence relation Eq. (B-22). Suppose {^J^) } 
is the minimal solution. Then, we have 



B 



(n) 
(1) 



(n-l) 



B 



(1) 



B 



(n) 
(1) 



fj{n)_ ^(n+l)_ ^(n+2)_ 
(n) ~(n-l)~(n) ~(n-2)~(n-l) 



y3(n)_ ^{n-l)_ ^(n-2). 



(1) 



-(l)-(2) ~(0)-(l) 



(B-24) 

(B-25) 
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From these equations, we can determine the ratios of aU the coefficients, B ^-j / B^-^ . 

Now, we determine the value of the two coefficients and B^^\ Because 

Eq. (B-19) must be satisfied for any value of x, we set x = 1 in Eq. (B-19) and 
obtain 



.B 



in) 



Im 



n=0 sBi^ 



{e)_fn + a 



1=0 s^im 



(B-26) 



Also, from the normalization condition Eq. (A- 7), we find 

1^ dxf^y ^ 4(-i)pK/3) 



P oo 
ni=0 



(«2) p{Q,/3) 
"2 



712=0 



(B-27) 



Because the Jacobi polynomials are orthogonal, we have 



1 



l + x 



13 



^ 2 r(n + a + l)r(n + /? + l)5n,,n2 
{2n + a + p+ l)r{n + l)r(n + a + /? + 1) 

Then, Eq. (B-27) reduces to 

2r{n + a + l)r{n + 13 + 



(B-28) 



oo 
n=0 



.A 



(n) 



.B 



(n) 



Im 



I] 



(2n + a + /3 + l)r(n + l)r(n + a + /3 + 1) 



(B-29) 



We can obtain the squares of s^^'^ and sB^!^'^ from Eqs. (B-26) and (B-29). Finally, 



we determine the signs of sA^m^ and sB\'^' by requiring that the sign of s5'°^(x) 
be the same as the sign of the spin-weight spherical harmonics in the limit au — > 0. 

In Table X, we list the eigenvalues of the spin-weighted spheroidal harmonics 
for q = ±0.9 from £ = 2 to 6 in the case Mu = 0.1. In Tables XI-XIII, we Ust the 
values of the spin-weighted spheroidal harmonics for cos 9 = —0.9, 0, 0.9, q = ±0.9 
and ^ = 2 - 6. 
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Table X. Eim for Mw = 0.1. 



I rn 




g = -0.9 






q = 0.9 






2 


2 


6 


2340859091 


X 


10" 


5 


7540002160 


X 


10° 


2 


1 


6 


1153623140 


X 


10° 


5 


8752937655 


X 


10° 


2 





5 


9957562220 


X 


10" 


5 


9957562220 


X 


10" 


2 


-1 


5 


8752937655 


X 


10" 


6 


1153623140 


X 


10" 


2 


-2 


5 


7540002160 


X 


10° 


6 


2340859091 


X 


10° 


3 


3 


1 


2175986227 


X 


-1 ,-.1 

10 


1 


1815913335 


X 


10 


3 


2 


1 


2116698236 


X 


10^ 


1 


1876700620 


X 


10^ 


3 


1 


1 


2057141841 


X 


10^ 


1 


1937158320 


X 


10^ 


3 





1 


1997300442 


X 


10 


1 


1997300442 


X 


10 


3 


-1 


1 


1937158320 


X 


10^ 


1 


2057141841 


X 


10^ 


3 


-2 


1 


1876700620 


X 


10^ 


1 


2116698236 


X 


10^ 


3 


-3 


1 


1815913335 


X 


10^ 


1 


2175986227 


X 


10^ 


4 


4 


2 


0141120497 


X 


10^ 


1 


9853070617 


X 


10^ 


4 


3 


2 


0105083654 


X 


10^ 


1 


9889090042 


X 


10^ 


4 


2 


2 


0069068104 


X 


10^ 


1 


9925093222 


X 


10^ 


4 


1 


2 


0033067556 


X 


10^ 


1 


9961086374 


X 


10^ 


4 





1 


9997075735 


X 


-1 ,-.1 

10 


1 


9997075735 


X 


10 


4 


-1 


1 


9961086374 


X 


10^ 


2 


0033067556 


X 


10^ 


4 


-2 


1 


9925093222 


X 


10^ 


2 


0069068104 


X 


10^ 


4 


-3 


1 


9889090042 


X 


-1 ,-.1 

10 


2 


0105083654 


X 


10 


4 


-4 


1 


9853070617 


X 


10 


2 


0141120497 


X 


10 


5 


5 


3 


0117824408 


X 


10^ 


2 


9877790649 


X 


10^ 


5 


4 


3 


0093445069 


X 


10^ 


2 


9901445845 


X 


10^ 


5 


3 


3 


0069155402 


X 


10^ 


2 


9925172191 


X 


10^ 


5 


2 


3 


0044953129 


X 


10^ 


2 


9948972039 


X 


10^ 


5 


1 


3 


0020835961 


X 


10^ 


2 


9972847731 


X 


10^ 


5 





2 


9996801598 


X 


10^ 


2 


9996801598 


X 


10^ 


5 


-1 


2 


9972847731 


X 


10^ 


3 


0020835961 


X 


10^ 


5 


-2 


2 


9948972039 


X 


10^ 


3 


0044953129 


X 


10^ 


5 


-3 


2 


9925172191 


X 


10^ 


3 


0069155402 


X 


10^ 


5 


-4 


2 


9901445845 


X 


10 


3 


0093445069 


X 


10 


5 


—5 


2 


9877790649 


X 


10^ 


3 


0117824408 


X 


10^ 


6 


6 


A 


2101144835 


X 




A 




X 




6 


5 


4 


2083478965 


X 


10^ 


4 


1912048117 


X 


IQi 


6 


4 


4 


2065910934 


X 


10^ 


4 


1928777542 


X 


IQi 


6 


3 


4 


2048439822 


X 


lOi 


4 


1945596361 


X 


10^ 


6 


2 


4 


2031064705 


X 


10^ 


4 


1962505533 


X 


10^ 


6 


1 


4 


2013784654 


X 


IQi 


4 


1979506008 


X 


IQi 


6 





4 


1996598734 


X 


lOi 


4 


1996598734 


X 


lOi 


6 


-1 


4 


1979506008 


X 


10^ 


4 


2013784654 


X 


10^ 


6 


-2 


4 


1962505533 


X 


10^ 


4 


2031064705 


X 


10^ 


6 


-3 


4 


1945596361 


X 


lOi 


4 


2048439822 


X 


10^ 


6 


-4 


4 


1928777542 


X 


10^ 


4 


2065910934 


X 


lOi 


6 


-5 


4 


1912048117 


X 


10^ 


4 


2083478965 


X 


10^ 


6 


-6 


4 


1895407126 


X 


10^ 


4 


2101144835 


X 


10^ 
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Table XI. Spin weighted spheroidal harmonics for Muj = 0.1 and cosO = —0.9. 



£ 


rn 


(J = -0.9 


q = 0.9 


2 


2 


4.3386981105 x 10"^ 


3.5950865567 x 10"^ 


2 


1 


-3.7087357243 x 10"^ 


-3.1984590567 x 10"^ 


2 





1.9411700165 X 10"^ 


1.7424257031 x 10"^ 


2 


-1 


-6.7728189805 x 10"^ 


-6.3276040875 x 10"^ 


2 


-2 


1.4469520014 x 10° 


1.4070158595 x 10° 


3 


3 


2.6579791433 x 10"^ 


2.3432796925 x 10"^ 


3 


2 


-2.3106765235 x 10"^ 


-2.0900105378 x 10"^ 


3 


1 


1.2376650604 x 10"^ 


1.1489471197 x 10"^ 


3 





-4.4858142218 x 10"^ 


-4.2770231715 x 10"^ 


3 


-1 


1.0509804998 x 10° 


1.0312667647 x 10° 


3 


-2 


-1.1695230417 x 10° 


-1.1936637613 x 10° 


3 


-3 


-9.1785502130 x 10"^ 


-8.8540840083 x 10"^ 


4 


4 


1.3964986285 x 10"^ 


1.2717039185 x 10"^ 


4 


3 


-1.2581385396 x 10"^ 


-1.1651738472 x 10"^ 


4 


2 


7.0862042637 x 10"^ 


6.6759828566 x 10"^ 


4 


1 


-2.7776127987 x 10"^ 


-2.6634215480 x 10"^ 


4 





7.5267562259 x 10"^ 


7.3543750538 x 10"^ 


4 


-1 


-1.2550656615 x 10° 


-1.2540767509 x 10° 


4 


-2 


6.8973736289 x 10"^ 


7.2640111281 X 10"^ 


4 


-3 


1.2548549058 x 10° 


1.2430573062 x 10° 


4 


-4 


4.9003837955 x 10"^ 


4.7271311939 x 10"^ 


5 


5 


6.8987908263 x 10"" 


6.4066680667 x 10"" 


5 


4 


-6.4686885018 x 10"^ 


-6.0782911246 x 10"^ 


5 


3 


3.8436104725 x 10"^ 


3.6551010190 x 10"^ 


5 


2 


-1.6262489225 x 10"^ 


-1.5656422649 x 10"^ 


5 


1 


4.9724643493 x 10"^ 


4.8496714304 x 10"^ 


5 





-1.0442277264 x lO" 


-1.0333677913 x 10° 


5 


-1 


1.2430710500 x 10° 


1.2553442781 x 10° 


5 


-2 


-9.9162402163 x 10"^ 


-1.3327302849 x 10"^ 


5 


-3 


-1.3250352372 x 10° 


-1.3300224425 x 10° 


5 


-4 


-8.7853110913 x 10"^ 


-8.6299962142 x 10"^ 


5 


-5 


-2.4419490099 x 10"^ 


-2.3603241859 x 10"^ 


D 


D 


O.zyo400l/OD X lU 


o.iuooiu/oyu X lu 


6 


5 


-3.2192050631 x 10"'^ 


-3.0549982369 x 10"^ 


6 


4 


2.0138768387 x 10"^ 


1.9279986867 x 10"^ 


6 


3 


-9.1318426706 x 10"^ 


-8.8215649196 x 10"^ 


6 


2 


3.0845665000 x 10"^ 


3.0079608194 x 10"^ 


6 


1 


-7.6021402310 x 10"^ 


-7.4896406387 x 10"^ 


6 





1.2535400677 x 10° 


1.2503499550 x 10° 


6 


-1 


-1.0079278456 x 10° 


-1.0282273588 x 10° 


6 


-2 


-4.8304664552 x 10"^ 


-4.5779972946 x 10"^ 


6 


-3 


1.1370699694 x 10° 


1.1530844521 x 10° 


6 


-4 


1.1913495835 x 10° 


1.1822538661 x 10° 


6 


-5 


5.2560056636 x 10"^ 


5.1488724800 x 10"^ 


6 


-6 


1.1736688454 x 10"^ 


1.1371714705 X 10"^ 
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Table XII. Spin weighted spheroidal harmonics for Muj = 0.1 and cos 9 = 0. 



£ 


rn 


g = -0.9 


q = 0.9 


2 


2 


4.1111211672 X 10"^ 


3.7950484443 x 10" 


i 


2 


1 


-8.0605525927 x 10"^ 


-7.7444328043 x 10" 


1 


2 





9.6770343965 x 10"^ 


9.6770343965 x 10" 


1 


2 


-1 


-7.7444328043 x 10"^ 


-8.0605525927 x 10" 


1 


2 


-2 


3.7950484443 x 10"^ 


4.1111211672 x 10" 


1 


3 


3 


5.8559507510 x 10"^ 


5.5982341302 x 10" 


1 


3 


2 


-9.3732214203 x 10"^ 


-9.3264155844 x 10" 


1 


3 


1 


7.1502176780 x 10~^ 


7.6307833481 x 10" 


1 


3 





3.8419145390 x 10"^ 


-3.8419145390 x 10" 


2 


3 


-1 


-7.6307833481 x 10"^ 


-7.1502176780 x 10" 


1 


3 


-2 


9.3264155844 x 10"^ 


9.3732214203 x 10" 


1 


3 


-3 


-5.5982341302 x 10"^ 


-5.8559507510 x 10" 


1 


4 


4 


7.1154532347 x 10"^ 


6.9134142570 x 10" 


1 


4 


3 


-9.8923696324 x 10"^ 


-9.9458839545 x 10" 


1 


4 


2 


5.0542053337 x 10"^ 


5.5505006446 x 10" 


1 


4 


1 


4.0134847826 x 10"^ 


3.4802911340 x 10" 


1 


4 





-8.3798585128 x 10"^ 


-8.3798585128 x 10" 


1 


4 


-1 


3.4802911340 x 10"^ 


4.0134847826 x 10" 


1 


4 


-2 


5.5505006446 x 10"^ 


5.0542053337 x 10" 


1 


4 


-3 


-9.9458839545 x 10"^ 


-9.8923696324 x 10" 


1 


4 


-4 


6.9134142570 x 10"^ 


7.1154532347 x 10" 


1 


5 


5 


8.1075684145 x 10^^ 


7.9470010304 x 10" 


1 


5 


4 


-1.0107000628 X lO" 


-1.0200385934 x 10° 


5 


3 


3.3710130936 x 10"^ 


3.8104293488 x 10" 


1 


5 


2 


6.0312055648 x 10"^ 


5.6888007421 x 10" 


1 


5 


1 


-7.6566919046 x 10"^ 


-7.8489902785 x 10" 


1 


5 





-2.5485570450 x 10"^ 


2.5485570450 x 10" 


2 


5 


-1 


7.8489902785 x 10"^ 


7.6566919046 x 10" 


1 


5 


-2 


-5.6888007421 x 10"^ 


-6.0312055648 x 10" 


1 


5 


-3 


-3.8104293488 x 10"^ 


-3.3710130936 x 10" 


1 


5 


-4 


1.0200385934 x 10° 


1.0107000628 X 10° 


5 


-5 


-7.9470010304 x 10"^ 


-8.1075684145 x 10" 


1 


D 


D 


o.yZ/ZyDolOU X lU 


8.7970600210 x 10" 


1 


6 


5 


-1.0179015837 x lO" 


-1.0287597521 x 10° 


6 


4 


1.9938800690 x 10"^ 


2.3714685226 x 10" 


1 


6 


3 


7.2764306367 x 10"^ 


7.0598308423 x 10" 


1 


6 


2 


-6.6414145769 x 10"^ 


-6.8995712931 x 10" 


1 


6 


1 


-2.7229278396 x 10"^ 


-2.3136795981 x 10" 


1 


6 





8.1612099685 x 10"^ 


8.1612099685 x 10" 


1 


6 


-1 


-2.3136795981 x 10"^ 


-2.7229278396 x 10" 


1 


6 


-2 


-6.8995712931 x 10"^ 


-6.6414145769 x 10" 


1 


6 


-3 


7.0598308423 x 10"^ 


7.2764306367 x 10" 


1 


6 


-4 


2.3714685226 x 10"^ 


1.9938800690 x 10" 


1 


6 


-5 


-1.0287597521 x 10° 


-1.0179015837 x 10° 


6 


-6 


8.7970600210 x 10"^ 


8.9272963150 x 10" 


1 
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Table XIII. Spin weighted spheroidal harmonics for Muj = 0.1 and cos^ = 0.9. 



£ 


rn 


q = -0.9 


q = 0.9 


2 


2 


1.4070158595 x lO" 


1.4469520014 x 10° 


2 


1 


-6.3276040875 x 10"^ 


-6.7728189805 x 10"^ 


2 





1.7424257031 x 10"^ 


1.9411700165 X 10"^ 


2 


-1 


-3.1984590567 x 10"^ 


-3.7087357243 x 10"^ 


2 


-2 


3.5950865567 x 10"^ 


4.3386981105 x 10"^ 


3 


3 


8.8540840083 x 10"^ 


9.1785502130 x 10"^ 


3 


2 


1.1936637613 x lO" 


1.1695230417 x 10° 


3 


1 


-1.0312667647 x 10" 


-1.0509804998 x 10° 


3 





4.2770231715 x 10"^ 


4.4858142218 x 10"^ 


3 


-1 


-1.1489471197 x 10"^ 


-1.2376650604 x 10"^ 


3 


-2 


2.0900105378 x 10"^ 


2.3106765235 x 10"^ 


3 


-3 


-2.3432796925 x 10"^ 


-2.6579791433 x 10"^ 


4 


4 


4.7271311939 x 10"^ 


4.9003837955 x 10"^ 


4 


3 


1.2430573062 x 10° 


1.2548549058 x 10° 


4 


2 


7.2640111281 X 10"^ 


6.8973736289 x 10"^ 


4 


1 


-1.2540767509 x 10° 


-1.2550656615 x 10° 


4 





7.3543750538 x 10"^ 


7.5267562259 x 10"^ 


4 


-1 


-2.6634215480 x 10"^ 


-2.7776127987 x 10"^ 


4 


-2 


6.6759828566 x 10"^ 


7.0862042637 x 10"^ 


4 


-3 


-1.1651738472 x 10"^ 


-1.2581385396 x 10"^ 


4 


-4 


1.2717039185 x 10"^ 


1.3964986285 x 10"^ 


5 


5 


2.3603241859 x 10^^ 


2.4419490099 x 10"^ 


5 


4 


8.6299962142 x 10"^ 


8.7853110913 x 10~^ 


5 


3 


1.3300224425 x 10° 


1.3250352372 x lO" 


5 


2 


1.3327302849 x 10"^ 


9.9162402163 x 10"^ 


5 


1 


-1.2553442781 x 10° 


-1.2430710500 x 10° 


5 





1.0333677913 x 10° 


1.0442277264 x 10° 


5 


-1 


-4.8496714304 x 10"^ 


-4.9724643493 x 10"^ 


5 


-2 


1.5656422649 x 10"^ 


1.6262489225 x 10"^ 


5 


-3 


-3.6551010190 X 10"^ 


-3.8436104725 x 10"^ 


5 


-4 


6.0782911246 x 10"^ 


6.4686885018 x 10"^ 


5 


-5 


-6.4066680667 x 10"" 


-6.8987908263 x 10"" 


D 


D 


i.lO/l/14/UO X iU 


l.i/ODDoo404 X lU 


6 


5 


5.1488724800 x 10"^ 


5.2560056636 x 10"^ 


6 


4 


1.1822538661 x 10° 


1.1913495835 x 10° 


6 


3 


1.1530844521 x 10° 


1.1370699694 x 10° 


6 


2 


-4.5779972946 x 10"^ 


-4.8304664552 x 10"^ 


6 


1 


-1.0282273588 x 10° 


-1.0079278456 x 10° 


6 





1.2503499550 x 10° 


1.2535400677 x 10° 


6 


-1 


-7.4896406387 x 10"^ 


-7.6021402310 x 10"^ 


6 


-2 


3.0079608194 x 10"^ 


3.0845665000 x 10"^ 


6 


-3 


-8.8215649196 x 10~^ 


-9.1318426706 x 10"^ 


6 


-4 


1.9279986867 x 10"^ 


2.0138768387 x 10"^ 


6 


-5 


-3.0549982369 x 10"^ 


-3.2192050631 x 10"^ 


6 


-6 


3.1033107890 x 10"'' 


3.2984531756 x 10"" 
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